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The spectral theorem of many-body Green's function theory relates thermodynamic correlations 
to Green's functions. More often than not, the matrix F governing the equations of motion has 
zero eigenvalues. In this case, the standard text-book approach requires both commutator and anti- 
commutator Green's functions to obtain equations for that part of the correlation which does not lie 
in the null space of the matrix. In this paper, we show that this procedure fails if the projector onto 
the null space is dependent on the momentum vector, k. We propose an alternative formulation 
of the theory in terms of the non-null space alone and we show that a solution is possible if one 
can find a momentum-independent projector onto some subspace of the non-null space. To do this, 
we enlist the aid of the singular value decomposition (SVD) of T in order to project out the null 
space, thus reducing the size of the matrix and eliminating the need for the anti-commutator Green's 
function. We extend our previous workm, dealing with a ferromagnetic Heisenberg monolayer and 
a momentum-independent projector onto the null space, to models where both multilayer films and 
a momentum- dependent projector are considered. We develop the numerical methods capable of 
handling these cases and offer a computational algorithm that should be applicable to any similar 
problem arising in Green's function theory. 

PACS numbers: 75.10.Jm, 75.70.Ak 



I. INTRODUCTION 

The spectral theorem in many-body Green's func- 
tion (GF) theory relates thermodynamic correlations to 
Green's functions, thus providing equations which, when 
iterated to self-consistency, allow the computation of 
the expectation values of the operators from which the 
Green's functions are constructed. Each Green's func- 
tion can be expressed in terms of an inhomogeneity of 
the equations of motion and higher-order Green's func- 
tions. Each higher-order GF can in turn be expressed 
in terms of yet another inhomogeneity and even higher- 
order Green's functions, and so on. Truncation of this 
infinite, exact hierarchy seldom occurs naturally. It is 
usually brought about through a decoupling approxima- 
tion, whereby the GFs of some order are approximated as 
linear combinations of lower-order functions which have 
already appeared in the hierarchy. This leads to a closed 
system of linear equations, the so-called equations of mo- 
tion, which relate the Green's functions to the inhomo- 
geneities. Anticipating results from the detailed exposi- 
tion in the next section, we may write the equation of 
motion in compact matrix notation: 

(wl-r)G = A. (1) 

Here, G is a vector whose components are the Green's 
functions, A is a vector of associated inhomogeneities. 
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and r is the matrix (in general unsymmetric) containing 
the coefficients obtained via the truncation. 

Each Green's function has poles at all the eigenvalues, 
uji, of the matrix T. As we shall show in detail in the 
overview section below, the spectral theorem associates 
a vector of correlation functions, C, to quantities at these 
poles, and thereby provides a route from the inhomogene- 
ity vector to the vector of correlation functions. In fact, 
this relationship can be expressed compactly|3| by a mul- 
tiplication of the vector A by a matrix constructed from 
the eigenvalues, LUi, and the right and left matrices of 
eigenvectors, R, L of the matrix F: 

C = RfLA, (2) 
LFR = O, (3) 

where is a diagonal matrix having the eigenvalues uji on 
the principal diagonal and £ is a related diagonal matrix 
with elements = (5y /(e''"'-l) , /3 being l/(fcT). Since 
C, A, and F all depend on the set of expectation values of 
the operators used in constructing the Green's functions, 
Eq. El can be solved by varying the expectation values 
until the left and right hand sides of the equation are 
equal. Usually, the equation holds in momentum space, 
so a Fourier transform to coordinate space (an integration 
over the momentum k) is also required; i.e. a set of 
integral equations has to be solved self-consistently. 

This procedure, while somewhat complicated to de- 
rive, is straightforward to apply, unless, as very often 
happens, the matrix F has a null space, i.e. there are 
eigenvalues of value zero. In that case, a naive use of 
Eq. 121 would involve a division by zero in the evaluation 
of £. The standard textbook procedure for handling the 
null space demands a knowledge of the anti-commutator 
GF in addition to the commutator GF from which Eq. [21 
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was derived (see references 0,0 and textbooks [1,0]). 
This leads to additional terms in Eq.|2 which restrict the 
GF method to the evaluation not of the complete vector 
C, but only to that part of C which lies in the non-null 
space of the matrix T. While this in itself is not neces- 
sarily a hindrance, it may become one if the projector 
onto the null space is dependent upon the momentum k. 
In this situation, even the modified equation Eq. |21 fails 
(see the following section); i.e. the standard textbook 
solution is in this case no solution. 

We propose here a new method that addresses this 
problem. It exploits the singular value decomposition 
(SVD) of the matrix governing the equations of motion 
in the many-body Green's function theory in order to 
treat the problems arising from the null space of the ma- 
trix. At each point in momentum space, the SVD effects 
a transformation to a smaller number of Green's func- 
tions having no associated null space, so that the corre- 
sponding correlations can be obtained from the spectral 
theorem in the straightforward manner of Eq. [5] This 
has the double advantage of both eliminating the need 
for the anti-commutator GF and reducing the size of the 
F-matrix, thereby also reducing the number of coupled 
equations needed to iterate to consistency. This idea was 
the subject of an earlier paper 0] on a ferromagnetic 
Heisenberg monolayer with single-ion anisotropy. That 
system is also solvable by the standard procedure because 
the projector onto the null space is momentum indepen- 
dent. As such, it is an inadequate model with which to 
demonstrate the effect of a momentum-dependent pro- 
jector. 

In the present paper, we treat as an example a model 
for Heisenberg thin films with exchange anisotropy, which 
does lead to a momentum-dependent projector. Note 
that the SVD of F is itself dependent upon k and so, 
while having the advantages mentioned above, does not 
solve the basic problem automatically. Nevertheless, it 
does allow one to find some SVD singular vectors (vide 
infra) which are in fact independent of k, thus provid- 
ing equations which can circumvent the problem of a 
momentum-dependent projector onto the null space. In 
connection with our method, there are a number of non- 
trivial numerical difficulties which arise and which be- 
come more acute when treating films with more than 
one layer. In order to confront these problems here, we 
extend our previous work 1] on the monolayer to multi- 
layer films and we describe a numerical procedure which 
surmounts these difficulties. 

In the following section, we give an overview of the es- 
sential equations of GF theory needed to explicate our 
new method, we state the precise nature of the problem, 
and we describe our proposal for solving it. The section 
after that outlines the numerical problems which arise 
and provides algorithms to solve them. Following that, 
we present the model for ferromagnetic Heisenberg thin 
films with exchange anisotropy as an illustrative exam- 
ple: we present an algebraic exposition of the model for 
a 3-layer film, which is easy to extend to any number of 



layers. A detailed study of the monolayer, for which all 
equations can be obtained analytically, reveals the struc- 
ture in the singular vectors from the SVD and the eigen- 
vectors of F, which harks back to the structure of the 
F-matrix. Some algebraic properties of the multilayers 
are then deduced from the structures found in the analy- 
sis of the monolayer. In the next section of the paper, we 
present some numerical results as an illustration. In the 
penultimate section we offer some remarks on how to use 
the new method to aid in the design of an efficient nu- 
merical algorithm. The last section contains a discussion 
and summary. 



II. OVERVIEW 
A. The standard formulation 

In this section, we supply all of the formulae neces- 
sary to make the paper self-contained and to allow the 
reader to understand and make use of the computational 
algorithm summarized in the discussion section. 

We start with the definition of the retarded Green's 
functions, which we shall use exclusively in this work: 

= {{Mt);B,{t'W, (4) 

where Q{t-t') is equal 1 for t > t' and for t < t' . Ai{t) 
and Bj{t') are operators in the Heisenberg picture, and 
i,j are lattice site indices. 77 — ±1 denotes anticommuta- 
tor or commutator GF's, (...) = Tr{...e->^'^)/Tr{e->^'^) is 
the thermodynamic expectation value, and Ti. the Hamil- 
tonian for the system under investigation. The GFs have 
a label a because for multidimensional problems, GFs 
are required for several different operators A and B. At 
this point, one need not be more specific. 

Taking the time derivative of equation and per- 
forming a Fourier transform to energy space, one obtains 
an exact equation of motion 

io{{A;B,))l^ ^ {[A,B,]X + {{[A,nUB,))l^. (5) 

Repeated application of the equation of motion to the 
higher-order Green's functions appearing on the right- 
hand side of Eq. \^ results in an infinite hierarchy of 
equations. In order to obtain a solvable closed set of 
equations, the hierarchy has to be terminated, usually 
by a decoupling procedure which, when restricted to the 
lowest-order functions, leads to the set of linear relations 

(([^„H]_;i?,))^^5]Fr^((A„;S,))^; , (6) 

/3,m 

where F is a matrix (in general unsymmetric) express- 
ing the higher-order Green's functions in terms of linear 
combinations of lower-order ones. 
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The lattice site indices can be eliminated by a Fourier 
transform to momentum space and the labels a can be 
suppressed by writing the equation of motion in compact 
matrix notation, where the vectors have components in- 
dexed by the labels a: 



{uj\ - r)G,, = A, 



(7) 



where 1 is the unit matrix, and the inhomogeneity vector 
has components A" = ([A, -B],,)". Note that G,, depends 
on energy and momentum (G^ = G^(aj,k)) and that 
A+i depends upon the momentum k, whereas A_i does 
not. 

It is now convenient to introduce the notation of the 
eigenvector method of reference 01 , since it is particularly 
suitable for the multi-dimensional problems in which 
many zero eigenvalues are likely to appear. One starts 
by diagonalizing the matrix F 



(8) 



where Jl is the diagonal matrix of N eigenvalues, 
lOt (t = 1, ...,N), No of which are zero and {N — No) 
are non-zero. The matrix R contains the right eigen- 
vectors as columns and its inverse L = R~^ contains 
the left eigenvectors as rows. L is constructed such that 
LR = RL = 1. Multiplying equation (TJ from the left by 
L, inserting 1 = RL, and defining new vectors Gri — LG^ 



and An 



LiA,j one finds 



(9) 



The crucial point is that each of the components r of this 
Green's function vector has but a single pole 



{Ari)T 



(10) 



This allows the standard spectral theorem, see e.g. [^|^, 
to be applied to each component of the Green's function 
vector separately: One can then define the correlation 
vector C = LCk, where Ck = (BA) is the vector of 
correlations associated with G,, (the index k is added to 
emphasize that one is in momentum space). 

For the commutator functions (ji = —1), the N — No 
components r for lOt 0, (C^)r (the upper index 1 refers 
to the non-null space) , are 



— lim 

27r -5^0 



doj 



1 



(5-i(w -I- id) - Q-i{uj - id))r 



(11) 



This equation cannot be used to define the Nq compo- 
nents of C° (the upper index refers to the null space) 
corresponding to w,- = because of the zero in the de- 
nominator. Instead, one must enlist the help of the anti- 
commutator Green 's functiondHllQ: 



(C°)ro = liin ^t^(^r;=+l)ro- 



The components of C°, indexed by tq, can be simplified 
by using the relation between the commutator and anti- 
commutator inhomogeneities, A+i = A_i -|- 2Ck. This 
yields, for ujra = 0, 



lim 



Ll>{A+i)-, 



1 



1)to 



2 (^^0 UJ — UJro 

i(L°(A_i + 2Ck)).o = (L"Ck)..„ 



(13) 



where we have used the regularity condition 0, L"A_i = 
0, which derives from the fact that the commutator 
Green's function is regular at the origin: 



lim L"(cj1 - r)G_i = L°A_i = 0. 



(14) 



Again, we use superscripts and 1 to denote the vec- 
tors belonging to zero and non-zero eigenvalues, respec- 
tively. The right and left eigenvectors and the correlation 
vectors may then be partitioned as 



R = (Ri R") L 



L" 



C 



(15) 



where the correlation vectors from equations 111|) 
and ini) are then = L^Ck and = f^L^A-i, and 
£i is the {N - No) x (TV - No) matrix with l/{e'^^- - 1) 
on the diagonal. 

Multiplying the correlation vector C from the left by 
R yields a compact matrix equation for the original cor- 
relation vector Ck: 



(l-R"LO)Ck = R^f^LiA. 



(16) 



In the above equation and for the rest of this paper, an 
inhomogeneity without a subscript rj always refers to the 
case rj = —1, i.e. A = A_i. Eq.^|is in momentum space 
and the coupled system of integral equations obtained 
by Fourier transformation to coordinate space has to be 
solved self-consistently. Usually one is interested only 
in the diagonal correlations; i.e. one has to perform an 
integration over k in the first Brillouin zone 

C = j rfkCk, (17) 

where the C without the subscript denotes the vector of 
diagonal correlations in coordinate space. 



B. The problem and a proposal for solving it. 



The quantity (1 — R'^L'^) functions as a projection op- 
erator onto the non-null space, and as such has no in- 
verse, so that Ck cannot be extracted from Eq. ^| This 
tells us that the Green's function method as formulated 
here can retrieve only part of the full correlation vector 
in coordinate space. 

In cases where the projector is momentum- 

(12) 

independent, which is the case in most of our previous 
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work 0, S Hi 01 it is not necessary to know 
the complete Ck, because one can take the projector 
outside the k-integration, solving the resulting equation 
self-consistently in coordinate space: 

(1 - R°L°)C = y dkRif^L^A. (18) 

The matrix elements of F, the inhomogeneity vector A 
and the correlation vector C in coordinate space depend 
only on the thermodynamic expectation values (magne- 
tizations and their moments), which are the variables in 
terms of which the system of equations H18|) is solved. 
This procedure was followed in most of our previous work 

Should, however the null-space projector depend on k, 
the above formulation cannot be used to solve for the 
expectation values because there is no expression for Ck 
available; i.e. the standard procedure fails. A way around 
this is to transform the Green's functions so as to elimi- 
nate the components lying in the null space. The spectral 
theorem then leads to a working equation in terms of a 
correlation vector which lies in the non-null space. 

The tool for finding the necessary transformation is 
the singular value decomposition (SVD) 0, ^| of the 
F-matrix: 

F = UWV. (19) 

The matrix W is a diagonal matrix whose elements are 
the singular values, which are > and U and V are 
orthonormal matrices containing the singular vectors: 
UU = 1 and VV = 1, where V denotes the transpose 
of V. Note that the F matrix is fully determined by the 
non-zero singular values and their corresponding singular 
vectors alone: 

F = UWV = UWV. (20) 

u and V arc N x (N — Nq) matrices obtained from U 
and V by omitting the columns corresponding to the zero 
singular values. The matrix w is the {N — Nq) x (N — Nq) 
diagonal matrix with the non-zero singular values on the 
diagonal. The remaining N x Nq matrices associated 
with the null space are denoted by Uq and Vq. The N x 
A'' matrices vv and VqVo are projectors onto the non- 
null space and null space of F, respectively. The sum 
of these two projectors spans the complete space: vv -f- 
vqVo = 1. It should be borne in mind that although 
R^L'^ also behaves as a projector onto the null space, 
one cannot identify Vq with R° or vq with L°, since R° 
and L'^ result from diagonalization of the unsymmetric 
matrix F. In fact we see from Eq. [201 that must lie 
in the space spanned by Uq so that (L'^u)wv — L°F = 
and similarly R° must lie in the space spanned by vq. 
Note that V and U are matrices of the eigenvectors of 
the symmetric matrices FF and FF, respectively; the 
eigenvalues of both of these matrices are the squares wf 
of the singular values of F. 



The crucial point is that the dimension of the equations 
of motion can be reduced by the number of zero singular 
values, which is equal to the number of zero eigenvalues 
of F, by applying the following transformations 

7 = vFv, (21) 

g = vG, (22) 

a = vA, (23) 

c - vCk. (24) 

This can be shown by multiplying Eq. ^ from the left 
by V and using the identity (recall that vv = 1) Fvv = 
uwvvv = F: 

v(tjl-Fvv)G = vA 
(t^l-vFv)vG = vA 

(c^l-7)g = a. (25) 

Here again, the eigenvector method can be used. Ma- 
trices 1 — L^v and r = vR^ diagonalize the 7-matrix, 
l7r ~ uj^ , where uj^ is identical to the matrix of non-zero 
eigenvalues of the full F-matrix (see Appendix A). The 
spectral theorem applied to the equation of motion with 
the matrix 7, which now has no zero eigenvalues, yields 
the equation for the correlations c = vCk in momentum 
space: 

c ^ r£ha, (26) 

where the {N — Nq) x {N — Nq) diagonal matrix £^ has 
the same elements as before, l/{e^'^^ ~ 1)- 

In order to determine the correlations in coordinate 
space, one has to perform a Fourier transform and then 
self-consistently solve the system of integral equations 

= y" dk(r£4vA - vCk) , (27) 

where A and Ck are the inhomogeneity and correlation 
vectors of the original problem, and the integration is 
over the first Brillouin zone. As pointed out in the intro- 
duction, Eq.|23still has the problem that the row- vectors 
in the matrix v are in general dependent on the momen- 
tum; however, for the model we have investigated, it is 
possible to find some row-vectors which are indeed inde- 
pendent of k. Only these rows in Eq. |23 may be used as 
working equations, for only then can those row-vectors 
in the last term on the rhs of Eq. |^ be taken outside 
the integral, allowing the term to be evaluated from the 
correlation vector in coordinate space: 

J dk VjCk vj y dk Ck = VjC, (28) 

where the index j labels one of the k-independent row- 
vectors. 

The above procedure formally solves the problem of 
the null space, but there remain two non-trivial problems 
born of the need to use numerical methods to obtain the 
vectors v and Vq in most cases. A precise statement 
of these problems and a description of the procedures 
needed to solve them are the subject of the next section. 
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III. METHODS FOR SOLUTION 

Here we describe in detail the numerical difficulties 
which arise because of the arbitrariness associated with 
any numerical determination of the singular vectors: vec- 
tors belonging to degenerate singular values are deter- 
mined only up to an orthogonal transformation and the 
phases of non-degenerate vectors are not fixed [T^. This 
arbitrariness may be removed by a smoothing procedure 
applied either to the SVD vectors as they are, (i.e. un- 
treated vectors) , or to vectors which have been previously 
subjected (optionally) to a labelling procedure. The 
somewhat lengthy description here will be shortened to 
a recipe given later in the discussion section. 



A. The numerical difficulties 

When the row vectors in v are determined numerically, 
they are unique only up to a sign change or, in the case 
of degenerate singular values, an orthogonal transforma- 
tion of the degenerate vectors: each time these vectors 
are computed anew, they are in effect rotated by a ran- 
dom amount with respect to vectors from the previous 
calculation; hence, even if the elements of the F-matrix 
are changed continuously (e.g. by varying the momentum 
k upon which they depend) , the computed vectors v from 
the singular value decomposition will not necessarily be 
continuous. This is also true of the vectors in matrices r 
and 1, which are obtained by a numerical diagonalization, 
but these occur in Eq. [SB] as factors separated only by a 
diagonal matrix having the same degeneracy structure 
as the eigenvectors, so that the arbitrariness in r and 1 
cancel each other. The vectors v, however, appear alone, 
so that vectors at neighbouring values of k will in gen- 
eral have arbitrary phases. If they differ by a change of 
sign, the integrands in Eq.[23exhibit discontinuities, pre- 
venting numerical evaluation of the integrals over k; we 
denote this problem by the term phase difficulty. If the 
vectors differ by an orthogonal transformation, the indi- 
vidual members of Eq. [23 will refer to different things at 
each value of k, and no meaningful system of equations 
can result; we shall call this the labelling difficulty, for 
reasons which will become apparent. Both of these diffi- 
culties can be overcome by rotating the vectors v among 
themselves to obtain a new set which spans the same 
space, is labelled, and has a smooth dependence on k. 
It should be emphasized that this rotation preserves the 
exact nature of the vectors v; they are not renormalized 
or adjusted in any other way whatsoever. 

At this point, it is pedagogically useful to consider a 
concrete example of the difficulties. To do this wc antici- 
pate the model of section IV for a 2-layer film, for which 
the matrix T has dimension 6 and the null space has 
dimension 2; hence, the two null vectors from the SVD 
exhibit both the labelling difficulty (because the vectors 
belong to degenerate singular values) and the phase diffi- 
culty. The first three components of these vectors refer to 



layer- 1 and the second three refer to layer 2. We calculate 
the dot product of each of the two vectors vq with the 
vector (1, 1, 1, 0, 0, 0) (lying fully in the space of layer-1) 
as a function of the momentum k along a line — ky 
in the first Brillouin zone. The numerical computation 
of the dot products suffers from both of the difficulties 
mentioned above, as shown in panel (a) of Fig. ^ where 
the full circles denote the dot product with the first vec- 
tor and the full triangles the dot product with the sec- 
ond vector. Here the words "first" and "second" refer 
to the order in which the vectors are returned from the 
subroutine calculating the singular value decomposition 
(the untreated vectors); they have no other significance. 
Clearly, the behaviour of the vectors as a function of 
k is unacceptable, for it would prevent our successfully 
evaluating the integrals numerically. Because of the way 
in which we calculate the integrals in Eq. [57| it is of 
the utmost importance to ensure that the vectors vary 
smoothly with k before the integral calculation is begun. 
In evaluating each integral, the range of k is divided into 
pieces and the contribution from each piece to the to- 
tal integral is summed. Then, the pieces are chopped 
into smaller pieces and the procedure repeated until the 
integral estimates (the sums) from two successive subdi- 
visions agree to within some error tolerance. This allows 
us to put the more effort into those regions where the in- 
tegrands are largest or change rapidly. This means that 
two successive integrand calculations may correspond to 
values of k far removed from each other. Thus, global 
smoothness is necessary; i.e. the phases of the vectors 
over the whole range of k must be fixed prior to the cal- 
culation of the integrals. This is achieved by a smoothing 
procedure, which we shall describe presently, but first we 
address the labelling problem. 



B. A labelling procedure 

It may be necessary to label the vectors v and vg either 
to ensure that vectors at neighbouring values of k refer 
to the same things or to construct vectors having some 
specific property (e.g. k-independence). It is instructive 
to note that, whereas the null space vectors Vq (as de- 
livered by the SVD) have no labels distinguishing them 
from each other, the vectors v are labelled by the singular 
values to which they belong, provided only that the lat- 
ter are not degenerate. Because degeneracies can arise, 
for some particular value of k say, we cannot rely on this 
labelling to protect us from the kind of discontinuities 
shown in panel (a) of Fig. ^ Now the singular values are 
of no importance in themselves: they serve only to sepa- 
rate the full space into a null space and non-null space. 
We are therefore free to take any linear combination of 
the vectors v among themselves or Vq among themselves. 
This freedom permits an assignment of unique labels of 
our own choosing. 

A very general approach is to generate block labels: 
the Green's functions are separated into labelled sets, so 
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FIG. 2: Full matrix V expressed as a sum of the diagonal 
reference matrix V^d and the off-diagonal blocks. The matrix 
of reference vectors also has a block structure and is the direct 
sum of the row-vectors belonging to each block. 



Momentum 



FIG. 1: Dot products of the two vo vectors of a bilayer 
(see the model described in section IV) onto the vector 
(1,1,1,0,0,0) lying in the space of layer-1 as a function of the 
momentum k. (a) untreated vectors, (b) vectors labelled with 
layer index, (c) vectors labelled and smoothed (full line); vec- 
tors smoothed but not labelled (dashed line). Calculations 
are for a 2-layer film with exchange energy J = 100, ex- 
change anisotropy D = 0.7, — 0.1, = 0, at temperature 
T = 90. (The reorientation temperature for the magnetiza- 
tion is Tr = 91). 



that the F-matrix is partitioned into blocks characterized 
by the row label index and column label index. Just 
how the blocks are chosen depends upon the model under 
consideration, e.g. in the model used to construct Fig.^ 
{vide infra), each block corresponds to Green's functions 
having the same layer index. A set of reference vectors, 
Vrcf , can then be constructed by finding the SVD of the 
associated matrix Fref which is just the blocked F-matrix 
with all off-diagonal blocks set to zero. This must be done 
so that the reference vectors also have a block structure: 
each vector has non-zero components only for the block 
to which it belongs, so that it may be labelled with that 



block index, as shown in Fig. El 

Fref = UrefW,efV,ef . (29) 

The set of reference vectors so constructed span the same 
space as all of the original untreated vectors (v and Vq 
together): VV ~ Vj-cfVi-of- It is convenient to define a 
block-label operator in terms of the reference vectors as 

Nb 

Fop:=Y.'^lS,mV% (30) 

where L{i) is some label for block z(for a specific choice, 

see Appendix B), Nb is the number of blocks, and V^*f 
is the set of reference vectors belonging to block i. The 
matrix of Pop in the basis of the singular vectors v (or, 
analogously, vq, if these are needed) of the full F matrix 
is 

P = vPopV = ^S,S„ (31) 

where Si := •\/-L(i)V^*fV This is equivalent to defining 
P as a product of overlap matrices expressed in terms of 
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a matrix of the reference vectors each multiphed by the 
square root of theh labels: 



P = SS, 



(1) 



(32) 



v/I(iVs)v(fJ-^ V. (33) 



See Fig. 121 for the meaning of the direct sum symbol ©. 
We now write S in terms of its singular value decompo- 
sition: 



S = LyZ, 
Zy^Z. 



SS 



(34) 
(35) 



i.e. Z diagonalizes SS, the matrix of the block-label op- 
erator. The singular values of S are the square roots of 
the eigenvalues of the block-label operator in the basis 
V. In other words, if each vector in v could be brought 
into coincidence with one of the reference vectors, then 
these vectors would have the same labels as the reference 
vectors. This does not, in general, occur, because the off- 
diagonal blocks of the F-matrix ensure that the reference 
vectors cannot all be singular vectors of the full matrix; 
hence, the squares of the singular values of S can be used 
as labels of transformed vectors, v^: 



VL = Zv. 



(36) 



The matrix of the block-label operator in the basis 
is , a diagonal matrix. This tells us that we have found 
a rotation (recall that ZZ = 1) that associates the vectors 
Vi as closely as possible with the reference vectors. If the 
block-label method is to be practicable, then the values 
y"^ should lie close to the labels L{i). 

In the concrete example of Fig. panel (5) shows the 
dot products again, this time computed from vectors Vo,l 
labelled with a layer-index. The erratic behaviour has 
now disappeared: it is clear that one of the null vectors, 
(the one whose dot product is nearly zero) is now asso- 
ciated largely with layer 2, while the other is associated 
with layer 1. Note that the discontinuities coming from 
the arbitrary signs of the vectors still remains. This can 
be cured with the smoothing procedure. 



C. The smoothing procedure 

At this point, it is convenient to drop the bold-face 
k, writing simply k instead, since we can, without loss 
of generality, discuss the methodology in terms of 1- 
dimensional integrals only. In fact, our models employ 
the simple cubic lattice, where the 2-D integrals can be 
reduced to 1-D integrals by exploiting the symmetry in 
the terms containing and ky. The first Brillouin zone 
is mapped onto the k-\me in the interval < fc < 1. Note 
that k is not simply the magnitude of k, but rather a pa- 
rameter which determines the value of all terms which 
depend on the momentum k. 

The smoothing procedure is akin to the labelling 
method in that untreated vectors V, which we shall here 



refer to as target vectors, are brought into as close a co- 
incidence as possible with the reference vectors, i.e. the 
target vectors are rotated to match the (fixed) reference 
vectors as well as possible. Here, however, the reference 
vectors are sets of standard vectors determined at var- 
ious reference values of k in the interval < fc < 1. 
As such, the reference vectors do not span exactly the 
same space as the set of target vectors at a neighbouring 
value of fc, so a slightly different procedure must be em- 
ployed. Note that in this section, the untreated vectors 
V (V = (v, Vq)) may or may not be labelled vectors V^, 
so we shall drop the subscript L for generality. If the 
vectors are not labelled, the method presented here must 
be powerful enough to cure the erratic behaviour shown 
in panel (a) of Fig. ^ not just the sign changes shown 
in panel (6) of the same figure. (See the dashed line in 
panel (c) of the figure for such a case.) 

We achieve the global smoothing in two stages. Prior 
to the integral calculation, we compute sets of well-placed 



reference vectors, {V 



(0., 

rcf ' ' 



1,. . .Nr} corresponding to 



points in the fc-interval. These allow us to generate by 
interpolation an appropriate reference vector anywhere 
on the fc-interval. Then, during the calculation of the in- 
tegrands in Eg. 12 71 we rotate the vectors v at a particular 
k so as to match as closely as possible the interpolated 
reference vectors appropriate for that fc. 

To construct the sets of reference vectors, we start by 
obtaining the vectors Vi-ef(fco) at some point ko, taking 
these to be primary reference vectors. Moving away from 
ko a small distance (e.g. a tenth of the fc— range), we cal- 
culate the vectors V(fci) at the trial point fci and from 
these, we construct a set of reference vectors Vrof (fci) by 
a procedure designed to find the best match of Vfcf (fci) 
with Vi.of(fco)- (An exposition of the methods to find a 
general rotation matching two vector spaces can be found 
in Appendix B.) The overlaps of corresponding vectors in 
the two sets are used as a criterion for accepting or re- 
jecting the trial point. If the point is accepted, points 
further away are tested until either some maximum al- 
lowed distance is reached, or a point fails the test. If the 
first trial point is rejected, the step-size is halved to get a 
new test point. In this way, a set of secondary reference 
points can be found linking the vectors over the whole 
fc-interval to the primary reference vectors. For the work 
reported here, sets of about 10 to 25 reference points 
sufficed to ensure a very large overlap (0.98, where 1.0 
is "perfect" ) of corresponding reference vectors at neigh- 
bouring fc-points. 

Once the sets of reference vectors are in place at the 
Nr reference points {ki,i = 0, 1, 2, . . . , A^,.}, smoothed 
vectors at any fc may be obtained in two steps: 1) a set 
of reference vectors at fc is obtained by interpolating be- 
tween the vectors at the reference points fc; and fc/i which 
bracket fc: fc/ < fc < fc/j. 2) The vectors v (untreated or 
labelled) are then rotated to be as close as possible to 
those in the interpolated reference set. 

The computation of the interpolated reference vectors 
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requires three steps: 

1. Weights are defined for the vectors at ki and kh'- 

2 TT / fc - fc; \ 

Wh = I -wi. 

2. A set of interpolated approximate reference vectors 

at k, Vref (fc) is obtained as a weighted sum of the 
vectors at ki and kh'- 

Vrof (fc) = W/Vrcf (fc/) + WhYrclikh). (37) 

3. The approximate vectors are orthonormahzed by a 
transformation found by diagonahzing their overlap 
matrix y: 

A = fyT, 

y-1/2 ^ TA^^/^t, 

V,ef = y-^^^^rci- 

The overlap matrix of the transformed vectors is clearly 
unity: 

VrefVref = TA-i/2(T3;T)A-i/2t, 
= TA-i/2(A)A-i/2t, 
= 1. 

The functional dependence of the weights chosen here 
ensures that the interpolated quantities are continuous 
and smooth at the reference points. We now have ref- 
erence vectors for the non-null space and the null space 

The procedure for matching the untreated (or labelled) 
vectors v with the interpolated reference vectors Vrcf dur- 
ing the integrand calculation is to seek a transformation 
Q that rotates the target vectors among themselves so as 
to achieve the best match, just as in the labelling proce- 
dure: 

vs = Qv. (38) 

Q is obtained via the singular value decomposition of the 
overlap matrix of the reference vectors with the target 
vectors: 

S := VrcfV, 

= /:xi, 
Q = cz, 

where x is the diagonal matrix of singular values of the 
overlap matrix, which are all very close to unity by con- 
struction. Q is indeed a rotation because 

QQ = Zi = 1, (39) 



the latter equality holding because the overlap matrix 
has no null space. The overlap of the transformed vectors 
with the reference vectors is almost the unit matrix: 

v,„fvs = SZC = CyiZZC = CxC « 1. (40) 

i.e. the new vectors are as close as possible to the ref- 
erence vectors. This procedure fixes the phases of the 
new vectors, because the transformation matrix Q is the 
product of the matrices C and Z which stem from a sin- 
gle SVD computation, so that the arbitrariness in the 
numerical determination of C and Z is lifted. 

To summarize: the untreated vectors can be both la- 
belled and smoothed by the transformation 

VLS = QiL = CZZv. (41) 

In our concrete example, the effect of the smoothing 
operation is seen in panel (c) of Fig. ^ The solid lines 
are the result of applying the smoothing to the labelled 
vectors in panel (b): Vo.Si = CZZ-Vq, where Vq are the 
untreated vectors labelled only by the singular values 0. 
The dashed lines result from smoothing the untreated 
vectors in panel (a): Vo,s = CZvq. 

In our example, we have not attempted to calibrate 
the primary vectors in any way, but this could be done if 
desired; however, this would still be a very crude way of 
labelling the vectors. A more important function of the 
primary vector is to ensure that all points in the region 
of a given point are related to one another, k is not 
the only parameter determining the elements of the F- 
matrix; the expectation values of the spin operators are 
others. In particular, if the derivatives of some quantity 
with respect to any of these parameters is needed, then 
it is a good idea to retain some standard vector with 
which to calibrate the primary vectors for all points in 
the neighbourhood of some given point. 



IV. HEISENBERG HAMILTONIAN WITH 
EXCHANGE ANISOTROPY 

A. Algebraic formulation 

In order to illustrate the labelling and smoothing 
procedures described here, we consider the model of 
reference|l4j. a Heisenberg Hamiltonian with an ex- 
change anisotropy for thin ferromagnetic films. In con- 
trast to the work in reference I'], the projector onto the 
null space here depends on the momentum. To show the 
importance of the labelling procedure, we need to extend 
the model of reference I'l to multilayer thin films. 

In what follows, a composite subscript of the type k,^ 
refers to the site k in layer k. We shall consider only near- 
est neighbour interactions in a simple cubic lattice struc- 
ture. The Hamiltonian is characterized by an exchange 
interaction with strength {Jk^ix > 0) between nearest 
neighbour lattice sites, a uniaxial exchange anisotropy in 
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the z-direction with strength {Dk^i^ > 0), and an ex- 
ternal magnetic field B = {B^ ,0, B^) confined to the 
reorientation plane of the magnetization: 



Eq. n]has the following structure: 



n= - 



2 Jk^iASk,S^^+ s^^s[-j 

<k,^lx> 
<k,^lx> 



(42) 



Here the notation — 



is introduced, where 
and l\ are lattice site indices and < kj^lx > indicates 
summation over nearest neighbours only. To fully de- 
scribe magnetization in the xz-plane for arbitrary spin 
5, Green's functions of the following type are required: 



GJf'X =« 5" ; {Sir{S7j" »> (43) 



where a is one of {+,— ,z} and m and n are integers 
which are determined by the spin S. The equation of 
motion for this Green's function is 



4 a,mn 



The generalized Tyablikov decoupling approximation [T^ 
is now applied to each Green's function on the rhs of the 
equation of motion, e.g. 



(45) 



Elimination of the site indices by means of a Fourier 
transformation to momentum space results in the ma- 
trix form of the equations of motion shown in Eq. ^ The 
components in the Green's function vector have the su- 
perscripts a and mn and two subscript layer-indices and 
depend on energy and momentum: G"^"™(aj, k). 



1. The 3-layer model 

We now specialize the exposition to the 3-layer film, 
since it is the smallest non-trivial example from which 
the extension to larger films is obvious. For this case, 



/En 

E21 

r32 



ojI 





E23 








V 









En E12 

E21 E22 

E32 













E23 



/Gn\ 
G21 
G31 
G12 
G22 
G32 

Gi3 

G23 
VG33/ 



Ell E12 
E21 E22 
E32 

/All\ 





A22 





VA33/ 



\ 








E23 



(46) 



Each of the entries in the 9x9 matrix T is in fact a 
3x3 matrix corresponding to a Green's function vec- 
tor with the same mn values and layer subscripts but 
with superscripts a = {+,— ,z} characterizing the vec- 
tor components. Using i as a layer- index, the diagonal 
matrix Ta is 





1 TT'J 

2^1 





1 TJX 

2^i 




(47) 



(44) where 



Hf 



Z, + {St)Mq~l^) , 

B'- + Duq{S^) + (J,,,+i + A:,,+l)(Sf+i) 

+ (J,;,,_l+A,.-l)(5f_i) , (48) 

B- + {Sf)Mq - 7k) + J^,^+l{St+,) + .h.^l{SU) , 

- {St)Duli. . 



For a square lattice and a lattice constant taken to be 
unity, 7k — 2(cosA:j; -I- cos ky), and g = 4 is the number 
of intra-layer nearest neighbours. The off-diagonal sub- 
matrices Tij for j — i ±1 are of the form 





■^Jij {Sf) 







(J.,+A,)(5f) 
-(J,,+A,)(Sf) 




(49) 

Unlike the elements of the F-matrix, the components of 
the vectors C and A are dependent on the values of m 
and n. For spin S — 1, the values of {m,n) required are 
(0, 1) and (1, 1). The values of the diagonal correlations 
Cii and An for each value of a are shown in Tabled The 
layer index i has been omitted for brevity. 

The block-diagonal structure of the matrix in Eq. 
allows us to write Eq. 1451 in terms of Green's functions 
and inhomogeneities labelled by a single layer index i. 
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TABLE I: Diagonal correlations and inhomogeneities for spin 
S = 1. 





a 


mn = (0, 1) 


mn — (1, 1) 




+ 




-2-(5"> + 3(5'"S") 


A 











z 


-{sn 






+ 






C 










z 


(5'") + (S-'S") 





Define singly-indexed Green's function vectors (here 3 x 
3 = 9 components), their corresponding correlations, and 
inhomogeneities for {i = 1, 2, 3}. 




The (9 X 9) matrices are independent of an index: 



r = 



rii 


ri2 





r2i 


r22 


r23 





r32 


r33 



(50) 



The big equations can now be replaced by 3 smaller equa- 
tions of motion and correlation equations i — 1,2,3: 

(c^i-r)G, = A, (51) 

(l-R°L°)Ci(k) = R^S^L^A,. (52) 

Applying the SVD to F yields six vectors v in the non- 
null space (each layer contributes one null vector). The 
corresponding six correlations c are obtained by multiply- 
ing the last equation by v and inserting vv + vqVo = 1: 



{0,efc,-efc}, where = ^/iFW^Th^^. The matri- 
ces of eigenvalues and eigenvectors of F are then 



LFR 



R 





1 



1 



(55) 



(efe - 7?^)#^ 2i?^i/^ 
-(efc -|-7?^)i7^ 2i?^i/^ 



Taking the first row of L from Eq. 1551 and the inhomo- 
geneity vectors for mn = {0,1} from Table ^ we find 
from the regularity condition L*^A = 0, 



IP 



(5^ 



(56) 



This implies that the ratio on the Ihs is not dependent 
upon the momentum. This can also be seen directly from 
the definitions in Eq. 0^ 



(5^) + Dq{S^ 



(57) 



From the same definitions it is clear that the ratio 
does depend upon momentum because differs from 
by a momentum-dependent term; hence the projector 
R^'L'^ onto the null space is momentum dependent. 

The singular vectors in U and V (see Eq. I20() may be 
obtained as the eigenvectors of the symmetric matrices 
rr and rr, respectively. The singular values in W are 
the square roots of the eigenvalues of these matrices. 



v(l-R"'L"')G,(k) =R^£:^L^(vv + voVo)^»- (53) 

Use of vR" = 0, L^vo = (see Appendix A), r = vR\ 
and 1 = L^v leads to 

va(k) =r£ilvA„ (54) 

which, after integration over k corresponds to Eq. 1271 

2. Analytical algebraic results for the monolayer 

We now investigate in detail the monolayer for spin 
5 = 1, for which many results can be obtained ana- 
lytically and reveal features which are pertinent to the 
structures found for the multilayers . 

The monolayer model leads to an equation of motion 
with the single diagonal block Fn and the vectors Gn, 
Cii, and All. For the remainder of this section, we 
shall drop the subscripts 11. The eigenvalues of T are 



W 

W22 



Wii 

M^22 





u = 



V = 



1 -H' 

72 ^W22 

1 -H' 
V2 V2W22 



V2W22 



_HZ_ \ 

2W22 * 
1W22 
W22 



/ -H' -H' 2g" \ 
/ yilVii V2W11 V2W11 

^ 



\ 



Wii 



~7P~ 

Wii 



(58) 

(59) 
(60) 

(61) 



(62) 



The matrix U has been blocked into u and Uq by the ver- 
tical line and V into v and Vq by the horizontal line. We 
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see here explicitly (by factoring out ff^ and applying the 
regularity condition) that the vector Lg, the first row- 
vector in the matrix L in Eq. 1551 is proportional to the 
vector Uq and is independent of momentum. In fact, the 
vectors in u arc also independent of momentum. Sim- 
ilarly, Ro is proportional to vq, but does depend upon 
momentum. As the momentum varies, vq also varies but 
remains in the plane containing the G^-axis and the line 
bisecting the axes and G~ . This implies that one 
of the vectors in the non-null space can be chosen to be 
perpendicular to this plane and therefore must be inde- 
pendent of momentum. This is the second row- vector in 
V in Eq. |B3 The first row- vector in V is also dependent 
upon momentum and lies in the same plane as the null 
vector. 

From these considerations, it is clear that, for the 
monolayer, only the second row- vector of V in Eg. 1621 is 
momentum-independent; hence, only the second row from 
Eq. (here specialized to the monolayer) can be used in 
the consistency equations. For spin S — 1, there are two 
expectation values which act as the variables which must 
be iterated to self-consistency in Eq.[23 (S^) and (S^S^). 
Since only one row of Ea . 1541 can be used, it is necessary 
to use this equation twice, each time with different values 
of mn (see TablcPl. The other correlations in Table|l]can 
be expressed in terms of {S^) and (S^S^) via the regu- 
larity condition L°A™" = 0, for m -I- n < 25' -I- 1 (for 
details see 0|)- 

In the multilayer case, each layer supplies an extra row 
in Eq. 1541 which can be used as a consistency equation, 
so that there are just enough equations to solve for the 
expectation values (Sf) and {SfSf) for each layer i. 

3. Algebraic properties of the multilayers 

While it is not practicable to perform the singular value 
decomposition algebraically for the multilayer films, it is 
nevertheless possible to deduce some properties of the 
singular vectors by examining the results for the mono- 
layer and combining these with the structure of the 3x3 
blocks of the matrix T in Eq. 021 

We start with the structure of the null vectors of the 
matrix F from Eq. |^ in a film of N layers. We may 
write a null row-vector in terms of the three components 
for each layer as 

vo = (4^\...,^r) (63) 

There arc N null vectors each with 3A^ components. We 
see from Eq. 1^3 that the components Vq^ of the null vec- 
tor for the ith layer considered in isolation have the form 
(ai,ai,bi). But this must also hold for each of the null 
vectors in Eq. ESI To see this, consider the F-matrix of 
Eq.[Sn| extended to N layers acting on a single null vector: 

Fvo = (64) 



This gives 3iV equations for the components of Vq, 3 for 
each layer. Denote the components belonging to layer i 
as (xi, UijZi) Now the structure of the 3x3 blocks in F 
is such that the first equations for each layer (i.e. equa- 
tions 1,4,7, . . .) involve only the components Xi and Zi 
for all z; i.e. the N equations suffice to determine the 
N ratios Xi/zi. The second equations for each layer (i.e. 
2,5,8,...) involve only yi and Zi and they are exactly 
the same as the first equations for each layer. The com- 
ponents Xi are therefore the same as the jji. But this 
implies, from the structure of the blocks in F, that the 
third equations for each layer are then automatically sat- 
isfied. The values of the Xi, yi, and Zi are then obtained 
from normalization. This means that each null vector, 
although perhaps dependent upon the momentum, must 
lie in the hyperplane in which the components belonging 
to layer i have the form (oi, a^, hi). This, in turn, means 
that the hyperplane in which all layer components have 
the form (c^,— Ci,0) must lie in the non-null space and 
be independent of momentum, since each vector of this 
form is obviously orthogonal to each null vector. 

Consider now special vectors of this form for which the 
components for layer i are (ci, — Ci,0) but all other com- 
ponents are zero, i.e. N of the 37V reference vectors of 
Eq. Viof . Clearly, these vectors are also in the non- 
null space and are orthogonal to each other. This is the 
reason why the labelling procedure picks out these par- 
ticular reference vectors as being identical with some of 
the vectors v. Of course, the rest of the vectors v are not 
in general identical to any of the remaining 2N reference 
vectors, because the null vectors do not in general have a 
layer-structure, implying that each of the remaining ref- 
erence vectors lies partially in the null space. All of these 
features are fully verified by the numerical calculations 
of the next section. 

In summary, we see that in the iV-layer film, there 
are only N components of the vector c from Eq. |^ (for 
a given value of mn) which are independent of momen- 
tum and also possess a layer-structure. The numerical 
labelling procedure is capable of delivering these vectors. 
This bit of serendipity is essential to the success of the 
entire method, since we only know how to determine the 
diagonal correlations, Ca, from the expectation values; 
hence, we require that the dot product of our vector with 
Ci (k) only involve the component Ca . This is ensured by 
the aforementioned property of those (labelled) vectors 
in V which are independent of momentum. In contrast, 
the momentum-dependent vectors in v do not share this 
property and therefore cannot be used in solving the con- 
sistency equations (|77|l . 



B. Numerical results 

We conclude our exposition with a few illustrative nu- 
merical calculations for the S = 1 Heisenberg Hamilto- 
nian model with exchange anisotropy of Eq. 1421 The 
exchange energy is J = 100, the exchange anisotropy is 
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TABLE II: 
film. 



Row- vectors V at = 0.0100113 for the 3-layer 
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1 






2 






3 
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+ 


- 


z 
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- 


z 


+ 


- 


z 
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FIG. 3: One component of a momentum-dependent vector v 
for the 3-layer film near the Reorientation temperature of the 
magnetization (T = 107) at the reference points ki. 



D = 0.7, and the field in the z-direction is zero. We show, 
as representative examples, some results of calculations 
for 1-, 3-, and 7-layer films. 

First, the labelled and smoothed row- vectors V (V = 
(v,Vo)) are shown in Table Hll for a 3-layer film with 
= 0.1, = at T = 107, which is very near the re- 
orientation temperature of the magnetization, Tr. Dis- 
played are the six vectors in the non-null space, v, and 
the three null vectors, Vq, at fc = 0.0100113. Vectors 
1,3, and 5 are independent of momentum and have a 
layer-structure: i.e. they may be identified with layers 
1,2, and 3, respectively. In contrast to this, vectors 2, 4, 
and 6 have contributions from each layer. Nevertheless, 
each of these vectors can be labelled by the layers 1,2, 
and 3, respectively, because the contributions from those 
layers are by far the greatest in each respective vector. 
The null-vectors are similar to vectors 2, 4, and 6: they 
vary with k, have no layer-structure, but may be labelled 
with a layer index. 

A component of one of the vectors v is shown in Fig.O 
at the reference points ki for one of the momentum- 
dependent vectors for the 3-layer film at a temperature 
T = 107 near the reorientation temperature of the mag- 
netization. Since there is rapid variation near fc = 0, 




0.02 0.04 0.06 0.08 0.1 
k 

FIG. 4: A single component of one of the vectors in v for the 
monolayer as a function of momentum k for increasing values 
of magnetic field B^: 5="° « < B==i . . . < ^ 13 The 
circles indicate the position of the reference vectors. 



the number of reference points there is denser than at 
larger k. At smaller temperatures, this behaviour is not 
as extreme. 

Fig. ^displays the most rapidly varying component of 
one of the momentum-dependent labelled vectors v in 
the non-null space as a function of the momentum k over 
the first tenth of the /c-range for the monolayer at T = 0; 
the curves remain roughly constant over the rest of the 
range. For increasing magnetic field in the x-direction, 
the curves become steeper near fc = 0, becoming 
nearly L-shaped for the largest field value (correspond- 
ing to {S^) — > 0). A greater density of reference vectors 
(positions are indicated by the full circles) is needed in 
this range. Nevertheless, the smoothing procedure is able 
to deliver acceptable vectors even under these extreme 
conditions. Note that all the vectors v in the non-null 
space are needed accurately, not just those that are in- 
dependent of momentum, because they are employed in 
reducing the size of the F-matrix by the transformation 
in Eqns.EIltoEl 

We now present the magnetizations calculated by solv- 
ing Eq.|23self-consistently as a function of magnetic field 
and temperature. Fig. O shows the magnetization (5^) 
(beginning at 1 at — 0), the magnetization (S^) in- 
creasing linearly from 0, the reorientation angle of the 
magnetization, 8, and the absolute value of the total 
magnetization for a spin S = 1 Heisenberg monolayer as 
functions of the applied magnetic field in the a;-direction, 
B^ at temperature T = 0. (Results at other tempera- 
tures are similar.) 

Magnetization components as a function of tempera- 
ture are shown for the 3-layer film in Fig. Magne- 
tizations for layers 1 and 3 are identical as required by 
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FIG. 5: Magnetization components (S^) (= 1 at = 0) 
and {S^} (solid lines), the reorientation angle of the magneti- 
zation, Q = (arctan |^^)/-| (dashed line), and the absolute 
value of the total magnetization (dash-dot chain) for a spin 
S=l Heisenberg monolayer as functions of the applied mag- 
netic field in the s-direction, _B^, at T = 0. 




60 80 1 00 

temperature 



140 



FIG. 7: Magnetization components in different layers (S'^)i 
and {S^)i (solid lines) and the reorientation angle of the mag- 
netization (dashed line), for a spin S = 1 Heisenberg 7-layer 
film as functions of the temperature, T, at =0.1 



A MORE EFFICIENT ALGORITHM 
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FIG. 6: Magnetization components in different layers {S^)i 
and {S^)i (solid lines) and the reorientation angle of the mag- 
netization (dashed line), for a spin 5* = 1 Heisenberg trilayer 
as functions of the temperature, T, at =0.1 



symmetry. Because the magnetic field component is 
so small, the 3 components in the a;-direction cannot be 
distinguished from one another on the scale of the dia- 
gram. 

Fig. demonstrates that it is also no problem to cal- 
culate the reorientation of the magnetization for a film 
consisting of 7 layers. Films with spin S > 1 can also be 
treated. In fact, all of our previous work on this mo del [T^ 
can be reproduced with the present method, so there is 
no need to present more examples here. 



In this section, we shall establish the connection be- 
tween the method developed here and the one we used in 
our previous work on ferromagnetic Heisenberg films with 
exchange anisotropy 14j. To do this, the working equa- 
tions are cast in a slightly different form that in some 
cases allows one to dispense with the smoothing proce- 
dure. This suggests that the new method serves not only 
as a tool to find appropriate consistency equations but 
also as a means of designing an efficient algorithm for 
solution. 

Our method has so far employed the matrix of sin- 
gular vectors v to transform the correlations in momen- 
tum space (c = vCk) so that c can be expressed via 
Eq. 1261 in terms of quantities referring only to the non- 
null space of the matrix T. Via the labelling procedure, 
the row-vectors v were then transformed so that a sub- 
set of them were independent of momentum; hence, the 
result in Eq. |5S1 could be exploited in the consistency 
equations Eq. [23 Let us now multiply the rhs of Eq. 
by the unit matrix vv; 



vCk = v{vr£hvA). 



(65) 



Consider now the consistency equation obtained from the 
momentum-independent vector Vj , which we take as the 
second row- vector V in Eq. (refer also to Eq. I28|l : 



dk Ck = v,C 



J dk Vj (vr£ilvA) 



(66) 



There is one such equation for each layer. Taken to- 
gether, these equations succinctly describe the method 
of solution in reference where, for each layer, a dif- 
ference of correlations (the Ihs) was related to the differ- 
ence between the first and second integrated components 



14 



of the rhs of the equation in round brackets. The idea to 
build this difference of components came from the em- 
pirical observation (by taking linear combinations of the 
equations of the standard procedure in the full space) 
that such difference terms make no contribution to the 
momentum-dependent term R°L°, thus eliminating the 
latter from the consistency equations. Here fEq. lt)6|) . tak- 
ing the difference is effected via multiplication by the vec- 
tors -Vj found systematically from the SVD coupled with 
the labelling procedure, thus replacing intuition with a 
systematic approach. 

The advantage of this way of solving the equations is 
that in this case there is no need to apply the smoothing 
procedure at each value of k, since the singular vectors 
V and V now both occur in the product in the integrand; 
i.e. an arbitrary sign arising from v is offset by that from 
V. We emphasize that smoothing would still be needed if 
there were degeneracies in the non-zero singular values, 
which was not the case in our previous work. Also, the 
labelling and smoothing procedures are required to find 
the appropriate in the diagnostic phase of the work. 

This suggests that the procedure in this paper can be 
used as a tool to find appropriate vectors vl5 via the 
labelling and smoothing procedures and then, once the 
vectors have been found, to design a more efRcient nu- 
merical procedure. We note here that the momentum- 
independent vectors for the exchange anisotropy model 
are particularly simple, being in fact just numbers. This 
need not always be the case: in general, these vectors 
will depend on the magnetizations and will vary as one 
moves through solution space (e.g. magnetizations as a 
function of temperature). In this case, a method of the 
type suggested above would require that the diagnostic 
tool be built into the numerical algorithm in order to find 
the appropriate vectors at each point in space. 



VI. DISCUSSION 

The main thrust of this paper is to demonstrate that 
the SVD is a very powerful tool for solving the problems 
which arise when the equations of motion matrix T has 
zero eigenvalues. In particular, the singular vectors in the 
non-null space, v, are used to reduce the size of the ma- 
trix that needs to be diagonalized and, at the same time, 
increase the numerical stability by eliminating the zero 
eigenvalues and the degeneracy associated with them, a 
fact of considerable practical importance, especially when 
some of the non-zero eigenvalues lie close to zero. In addi- 
tion, the suitably labelled and smoothed singular vectors 
serve to transform the consistency equations so as 
to deliver the correlation functions in coordinate space 
required for the solution of these equations. 

The method outlined here can be viewed as a recipe for 
direct numerical calculation . The recipe is quite simple: 

1. At the current point in solution space, generate a 
sufficient number of reference vectors Vi.cf(fci) at 



suitably chosen knots ki. 

2. At each desired value of fc, 

(a) Obtain the untreated singular vectors v in the 
non-null space from the SVD of the F-matrix. 

(b) If necessary, find labelled vectors from these: 
vl = Zv. If no labelling is desired, set Z = 1 
at this point. 

(c) Find the set of smoothed reference vectors 
appropriate to k and carry out the smoothing 
transformation: v^s — CZTiv. 

(d) Use the smoothed vectors to reduce the 
F-matrix by eliminating the null space. 
Diagonalize the reduced matrix to get the 
eigenvalues and eigenvectors: l7r = 

(e) Select the appropriate (i.e. properly labelled 
momentum-independent) vectors v to gen- 
erate the set of consistency integrands of 
Eq.EZI r£lvA- vCk 

3. Finally, solve the consistency equations Ea.lTTIfe.g. 
with the method described in Q). 

This recipe cannot be so general as to apply to any 
problem, for the labelling is dictated by the properties 
of the F-matrix, which will depend upon the nature of 
the physical model and the approximations used in the 
decoupling leading to the F-matrix. The method can also 
fail if it is not possible to find momentum-independent 
vectors. 

The method may also be useful as a diagnostic tool 
to aid in the design of an efficient computation method. 
Once the consistency equations have been found by our 
systematic procedure. Eg. 1661 can be used as a more effi- 
cient working equation, since it usually does not require 
the smoothing procedure. In other words, our procedure, 
used diagnostically, serves to identify the allowed consis- 
tency equations. We are then free to manipulate these to 
form a new set of allowed equations having simpler nu- 
merical requirements, thus leading to an optimal solution 
of any particular solvable problem. 

APPENDIX A: REDUCTION OF THE F-MATRIX 
WITH THE SINGULAR VALUE 
DECOMPOSITION (SVD) 

In Eq. 1211 the matrix F is reduced in dimension by 
the number of zero eigenvalues, Nq, via a transformation 
by the singular vectors spanning the non-null space, v. 
It is very important that the eigenvalues of the reduced 
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matrix, 7, be identical with the non-zero eigenvalues of 
r, oj^. In this appendix, we show under which conditions 
we can expect this to be the case. 

While the null eigenvectors of F, L'' and R°, lie com- 
pletely in the spaces spanned by UqUo and vqVoj respec- 
tively, the vectors belonging to the non-zero eigenvalues, 
and R^, may have components in both the null space 
and non-null space, (e.g. vR° — but vqR^ ^ in 
general). In what follows, we assume that the left and 
right eigenvectors of F have been constructed to be or- 
thonormal: LR = 1. 

Because L and R diagonalize F, we have, inserting 
F = Fvv and vv -I- voVq = 1, 

n = LFR, 
(o^^O = (li)(vv + voVo)Fvv(R°,R1), 

= ( ) , (Al) 

V e21 622 y ^ ^ 

where the individual matrix blocks are 

en = L"(vv + voVo)FvvR°, 

= L"FR", 

ei2 = L"FvvR\ 

= L"FR\ 

621 = L^(vv + voVo)rvvR°, 
= L^FR", 

622 = Li(vv + voVo)rvvR\ (A2) 

Since FR" = and L°F = 0, the blocks en, ei2, and e2i 
are zero as required. The other block is the sum of two 
terms: 

622 = LVvFvvR^ + L VoVoFR\ 

= l7r + LivoVoFRi, (A3) 

where we have defined 1 L^v and r — vli^. If the 
second term in Eg. IA3I were zero, then the vectors 1 and 
r would diagonalize 7 yielding the eigenvalues as de- 
sired. This condition is fulfilled if the eigenvectors R° 
span the null space, for we can then express the vectors 
Vq in terms of the R°, each of which is orthogonal to 
by construction (L^R'' = 0). In all of our calculations, 
this was indeed the case and was checked numerically at 
every reduction of F. 

APPENDIX B: MATCHING OF TWO VECTOR 
SPACES 

In the context of this paper, we are usually concerned 
with rotating a set of target vectors (which we do not 
want to disturb in any other manner) so as to match 
as closely as possible a set of (approximate) reference 
vectors which serve as labels or calibration vectors. The 



method outlined below treats this problem in a general 
way and shows how the SVD can effect this matching in 
an optimal way. Note that the notation in this appendix 
is independent from that in the main body of the paper. 

It often happens that there are two sets of A^- 
dimensional vectors which are related by a general ro- 
tation (i.e. a rotation about an arbitrary axis) and that 
we wish to find this rotation. If each of the two sets are 
subsets of the whole space, it may be that each set spans 
a slightly different subspace. This is the problem that 
we most often meet here: the two subspaces do not over- 
lap exactly but have a high degree of overlap. In this 
case, we wish to find the "best" general rotation in the 
sense that each vector of one set is brought into as close 
a coincidence as possible with some vector of the other 
set. 

Let r — 1 , . . . , Nii label the Nr column vectors of 
length N , the reference vectors, in the matrix R of di- 
mension N X Nfj. Similarly, let i = 1, . . . , Nt label the 
Nt column vectors of length N, the target vectors, in the 
matrix T of dimension N x Nt ■ We assume here that the 
vectors in R and T are orthonormal among themselves 
(RR — 1 and TT = 1.) We wish to find a rotation 
which, when applied to the target vectors, brings them 
into closest coincidence with the reference vectors. In 
order to ensure that there are no systematic degenera- 
cies in the singular values that we will calculate, we first 
associate the vector R, i (i.e. the column i of the ma- 
trix R) with the unique label L{i) = Nji + 1 — i (i.e. 
a label which decreases as i increases. Then, we define 
a weighted overlap matrix whose Nji x Nt elements are 
given by 

Srt ■■= /CWRr.T.t. (Bl) 

The SVD of the matrix S is then 

(B2) 

The matrix of the label operator in the basis of target 
vectors is 

'PntxNt = SS. (B3) 

Writing S in terms of the SVD in this equation shows 
that the matrix Z diagonalizes P to produce eigenvalues 
y^, so that the required rotation is just Z. i.e. 

T' = TZ. (B4) 

The justification for this transformation is that the ma- 
trix of the label operator in this new basis is diagonal: 

P' = S'S'=y2. (B5) 

If Nif — Nt and the reference and target vectors lie 
in the same s ubspa ce, then the singular values yi will be 
identical to L(i) provided that the singular values are 
numbered in decreasing order, the usual convention; this 
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is the reason for the particular choice of labels L{i) above. 
If the reference and target vectors span slightly different 
spaces, then there will be some discrepancy between the 
singular values and the square roots of the labels. 

If the number of reference vectors is less than the num- 
ber of target vectors, then some of the singular values 
of S will be zero. This means that a subset of the ro- 
tated target vectors will be as closely coincident with the 



(smaller) number of reference vectors as possible and the 
rest of the (rotated) target vectors will be arbitrary and 
span the (null) space of S left over. The arbitrariness of 
the vectors is due to the degeneracy of the null space. 

If there are more reference vectors than target vectors, 
then each (rotated) target vector will be well-defined but 
the set of target vectors will not span the space of the 
reference vectors. 
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